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We use the quantum Monte Carlo (QMC) techniques to calculate the static structure function 
S(q) of a one-component ion lattice and use it to calculate the thermal conductivity n of high- 
density solid matter expected in the neutron star crust. By making detailed comparisons with the 
results for the thermal conductivity obtained using standard techniques based on the one-phonon 
approximation (OPA) valid at low temperature, and the multi-phonon harmonic approximation 
expected to be valid over a wide range of temperatures, we asses the temperature regime where S(q ) 
from QMC can be used directly to calculate k. We also compare the QMC results to those obtained 
using classical Monte Carlo to quantitatively asses the magnitude of the quantum corrections. We 
find that quantum effects became relevant for the calculation of k at temperature T < 0.3 Lip, where 
flp is the ion plasma frequency. At T ~ 0.1 Up the quantum effects suppress k by about 30%. The 
comparison with the results of the OPA indicates that dynamical information beyond the static 
structure is needed when T < 0.1 flp. These quantitative comparisons help to establish QMC as a 
viable technique to calculate k at moderate temperatures in the range T = 0.1 — 1 flp of relevance 
to the study of accreting neutron stars. This finding is especially important because QMC is the 
only viable technique so far for calculating re in multi-component systems at low-temperatures. 

PACS numbers: 97.60.Jd, 26.60.Gj, 95.30.Cq 


I. INTRODUCTION 

Observations of transient phenomena in accreting neu¬ 
trons stars and magnetars BI2 have motivated recent 
attempts to model the thermal evolution of the outer¬ 
most regions of the neutron star called the crust mi- 
These studies have shown that the thermal conductivity 
of the crust plays a very important role in shaping tem¬ 
poral structure of x-ray emission from these systems. In 
this study we revisit the calculation of the thermal con¬ 
ductivity of the outer crust of the neutron star, where the 
typical densities are in the range of 10 s —10 11 g cm -3 and 
temperatures are expected to be in the range of 10' — 10 9 
K. Under these conditions nuclei are pressure ionized and 
form a solid phase at low temperature. In contrast, elec¬ 
trons are relativistic, weakly coupled and very degener¬ 
ate. For these reasons electrons dominate the thermal 
conductivity of the neutron star crust. 

Electron conduction in the neutron star crust is lim¬ 
ited by electron-ion scattering. Ions in the crust have 
large atomic number (25 < Z < 40), and ionic structure 
and dynamics are strongly correlated by Coulomb inter¬ 
actions. Consequently, the amplitude for electron scat¬ 
tering off different ions interfere. Accounting for such in¬ 
terference under arbitrary ambient conditions in a multi- 
component plasma (MCP) is a challenging many-body 
problem. At low temperature and for the case when the 
ground state is a simple solid containing only one ion 
species, pioneering works by Flowers and Itoh {7 and 
by Yakovlev and Urpin [5] have shown that the one- 
phonon approximation (OPA) adequately describes the 
collective response of the ions which is needed to cal¬ 
culate the electron scattering rate. More recent studies 
have greatly improved the description of electron scatter¬ 


ing in a one component plasma (OCP) by including the 
effects of the multi-phonon processes with the harmonic 
approximation u m- It is now possible to calculate 
the electron thermal conductivity of an OCP over the 
full range of temperatures of interest to neutron star as¬ 
trophysics. For MCP, classical molecular dynamics has 
proven to be useful and is expected to provide a quanti¬ 
tive description of the thermal conductivity in the high 
temperature limit when T > Qp mi- However, neither of 
these methods are suitable to describe electron scattering 
in multi-component systems at low temperature. This 
observation is the primary motivation for our study here 
where we establish that Quantum Monte Carlo (QMC) 
methods, which are simply generalizable to MCP, can be 
used to calculate the thermal conductivity at the rela¬ 
tively low temperatures of interest in astrophysics. 

Quantum Monte Carlo (QMC) methods have been suc¬ 
cessful in determining the ground state properties of di¬ 
verse strongly correlated many-particle systems and can 
be adapted to study MCPs at low temperature when 
quantum effects become relevant. In this study, which 
is a first step towards developing a QMC approach to 
MCPs, we use the QMC method to study the ion-ion 
correlations in an OCP at low temperature and calcu¬ 
late the thermal conductivity. We also use the dynamic 
matrix to calculate the detailed phonon spectrum at low 
temperature which is in turn used to calculate re. These 
calculations allow us to make quantitative comparisons 
between the results obtained with the Monte Carlo tech¬ 
niques and those obtained with various approximations 
used in the literature. These comparisons, allow us to re¬ 
assess various theoretical approaches in some detail and 
provide useful new insights about the role of quantum 
effects, the low-energy strength of the ion response func- 
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tion, Bragg scattering and multi-phonon effects. 

The rest of the paper is organized as follows. In Sec¬ 
tion [TT] we review the basic properties of an OCP and 
define the relationship between the electron conductiv¬ 
ity and the ion dynamical structure function S(u>, q). In 
Section |III| we discuss the phonon spect rum and use it 
to calculate S(u;, q) in OPA. In Section IV we describe 
and use both classical and quantum Monte Carlo meth¬ 
ods to calculate the static structure function S(q) and 
compare these results to those obtained in OPA. In Sec¬ 
tion |V] we calculate k using different approaches and ap¬ 
proximations, and we identify the temperature regimes 
where these approximations are valid. We also outline 
the strategies to extend QMC calculations to MCPs in 
future work. Throughout this paper we adopt the natural 
physical units with h = c = k# = 1- 


II. THERMAL CONDUCTIVITY OF THE 
OUTER CRUST 

In the simplest scenario the outer crust at given den¬ 
sity is composed of cold catalyzed matter of a single ion 
species. The mass number A and charge Z of the ions 
are density dependent and are determined by minimiz¬ 
ing the total energy of the system. The ground state of 
such matter is a strongly correlated OCP with bare nu¬ 
clei immersed in a degenerate and weakly coupled elec¬ 
tron gas. The characteristic energy of the electron is set 
by its Fermi momentum 

Pf = (37r 2 n e ) 1/3 

where n e is the number density of electrons, and pro is 
the mass density in units of 10 10 g cm -3 . For the den¬ 
sities of interest to our study, which are typically in the 
range of 10 8 — 10 11 g cm -3 , pp is much larger than elec¬ 
tron mass m e , and it is a good approximation to treat 
electrons as ultra-relativistic. In contrast, ions are heavy 
and correlated. One characteristic energy is set by the 
ion plasma frequency 

f AnZ 2 e 2 ni \ 1/2 

np = V M ) 

M 2 .9X10» K )(|) (£)->, (2 , 

where ni = n e /Z is the ion density, M ss Am p is the mass 
of the ion with m p being the proton mass, and e 2 ~ 1/137 
is the fine structure constant in natural units. 

The typical Coulomb energy is of order if 2 e 2 /a, where 

“ = (^) “ <147fm) (si) f '«' 3 (3) 


is the inter-ion distance. The temperature T provides a 
measure of the “extractable” kinetic energy or the ther¬ 
mal energy of the ions. The ratio between the Coulomb 
energy and the thermal energy of the ions is a measure of 
the importance of interactions in the plasma and defines 
the dimensionless Coulomb parameter 


r = 


Z 2 e 2 
aT ' 


( 4 ) 


In a weakly coupled plasma T <C 1, and Coulomb inter¬ 
action can be studied within perturbation theory. Nu¬ 
merical simulations of the OCP have shown that ions 
crystallize into a solid state when F > T m « 175 E2HI2I- 
The melting temperature of the solid 


Tyr 


Z 2 e 2 

aT m 




2/3 


- 1/6 

PlO 
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can be correspondingly defined. We note that electron 
screening modifies the Coulomb potential generated by 
the ion at large distances. The modified potential is 


V{r) 


^ e p -rk T F 

r 


( 6 ) 


where 




( 7 ) 


is the Thomas-Fermi (screening) momentum. Because 
k^pa < 1 for the densities of interest, screening will not 
greatly alter the nearest neighbor ion-ion interaction, and 
the Coulomb parameter of the OCP without screening 
continues to provide a reasonable measure of the strength 
of interactions and the melting temperature. In this pa¬ 
per we have chosen to present results at fiducial densities 
10 10 and 10 11 g cm -3 , labelled as LD and HD, respec¬ 
tively. The chemical compositions in these two cases are 
chosen according to Page and Reddy [5] for catalyzed 
matter in the outer crust. We list the physical condi¬ 
tions of the two cases in Table [i] 

The dynamic structure function in the crystalline state 
is given by 


S(w,q) = ^ E e_iq ' (Ri_Rj ) 

i,j =1 

X (8) 

where N is the total number of ions in the crystal, R; is 
the equilibrium position of the ion on the *’th site, u.;(f) is 
the displacement of this ion at time t, and (■ ■ - )t denotes 
the thermal average. 
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TABLE I. Key parameters for cold catalyzed matter in neutron star crust at two fiducial densities. The compositions are 
chosen according to Page and Reddy 


Name 

Lattice 

P (g cm 3 ) 

Ion 

a (fm) 

Q P (10 9 K) 

T m (10 9 K) 

Pf 1 (fm) 

k^p (fm) 

fco 1 (fm) 

LD 

BCC 

To™ 

34 Se 

149 

0.32 

0.75 

24 

249 

62 

HD 

BCC 

10 11 

igNi 

68 

0.87 

1.11 

12 

121 

28 


The electron thermal conductivity k is mainly limited 
by electron-ion scattering and can be written as 


n 2 Tn e 
3e F ^ K ’ 


( 9 ) 


where e F ~ p F is the electron Fermi energy, and effective 
collision rate is [7j 

o rZpF 

Uk= 3 Z(2 tt) 3 J 0 dqq3 \ v ^ 2S '^' ( 10 ) 

In the above expression 

\ v (q)\ 2 = e 2 \V(q)\ 2 (l-^j (11) 


is the square of the scattering matrix element for 
electron-ion interaction with momentum exchange q, and 
the screened Coulomb potential generated by the ion in 
momentum space is 


where 


V{q) 


AirZe 
e(q) q 2 ’ 


*(*) = ! + k ~f 


( 12 ) 


(13) 


is the static dielectric function in the Thomas-Fermi ap¬ 
proximation. The effects due to ion-ion correlations on 
electron scattering are included in Eq. (101 through 


S' K (q)= dw (S'(uhq))q w K (u/T,q), (14) 


where S'(to, q) is the inelastic part of the dynamic struc¬ 
ture function, (• • • )^ is the average over the direction of 
unit vector q = q /q, and 


w K (z = u/T,q) 



(15) 


We note that the elastic Bragg scattering does not con¬ 
tribute to the conductivity because it has been accounted 
for in the ground state which leads to the electronic band 
structure. 

Because the response at high frequency w flp can¬ 
not involve collective motion of the ions, we expect that 
most of the strength of the dynamic response will reside 
at energies that are comparable to flp. Therefore, when 


T > Up it is a good approximation to retain only the 
leading terms of w K (z, q) in z = ui/T in the integrand in 
Eq. (14), and the static approximation S' K (q) « S'(q) is 
valid, where 



dw (S'(o;,q))-, 


(16) 


is the inelastic part of the static structure function. At 
very low temperature T <C Up, however, the exponen¬ 
tial factor l/(e - “/ T — 1) in w K (z,q) dominates. In this 
limit the static approximation breaks down, and S' K (q ) <C 
S'(q). Between these two temperature limits two com¬ 
peting factors in w K (z,q) dominate in different ranges of 
q. For large-angle scattering (with large q values) the ex¬ 
ponential factor still dominates, and S' K (q) < S'(q). But 
for small-angle scattering (with small q values) the factor 
Pp/q 2 can dominate, and S' K (q) > S'(q). We will deter¬ 
mine these temperature limits for OCP using OPA in the 
next section. 


III. PHONON SPECTRUM AND THE 
DYNAMICAL RESPONSE 


At low temperature the characteristic distance scale 
for ion motion in lattice is 

Huk) 1 ' 2 (i) W2 ^ < 17 > 

which is much shorter than the inter-ion distance a. Un¬ 
der these conditions the restoring force on the ion is 
quadratic in the displacement Uj(i), and the detailed 
phonon spectrum can be calculated by using the dynamic 
matrix mm 


N 

D(k) = 2 sin 2 

4 = 1 


(k-R t \ 

'd 2 U(x)' 

V 2 ) 

<9x<9x T 


(18) 


The phonon frequencies w s (k) (s = 1,2,3) are obtained 
by solving the eigenvalue equation Mw 2 (k) — D(k) = 0, 
and the corresponding normalized eigenvectors e s (k) are 
the phonon polarization vectors. In the long wavelength 
limit, phonons have a linear dispersion relation 

w 8 (k) =c s (k)k + 0(k 2 ), (19) 


where c s (k) is the sound speed of the phonon mode in 
propagation direction k = k /k. Generally speaking, po¬ 
larization vectors e s (k) are neither parallel nor perpen¬ 
dicular to k. However, for ka < 1, two of the phonon 
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FIG. 1. The dispersion relations of the photon modes of the Coulomb lattice. Left panel: The phonon frequencies u> as func tion s 
of wave number k. The dashed and solid curves are for the transverse and longitudinal modes described by Eqs. (201 and (221, 
respectively. The shaded region gives the range of the eigenvalues of the dynamic matrix in Eq. (181. Right panel: The shaded 
region gives the range of |e s ■ k| for the eigenvectors e s (k) of the dynamic matrix. The results in both panels are computed for 
the LD matter. 


modes in a cubic lattice are approximately transverse, 
and the third mode is approximately longitudinal. 

In Fig. |I] we show the phonon dispersion relations and 
polarization of a body centered cubic (BCC) lattice which 
are calculated from the dynamic matrix. The lower fre¬ 
quency modes in the left-panel of Fig. |T] correspond to 
the modes that are mostly transverse with e s • k ss 0, 
and the higher frequency modes are mostly longitudinal 
with |e s • k| ss 1 . 

Despite the relatively large spread of velocities associ¬ 
ated with the transverse modes it is often useful to rep¬ 
resent them by an “average velocity” denoted as C(. This 
is typically defined by taking the limit of k —> 0 in which 
case 


c t 


fco 


0.0031 



-2/3 


1/6 
P 10 > 


( 20 ) 


where a ~ 0.39 according to Chabrier et al. and 


At low temperature it is useful to write the dynamic 
structure function as a sum of the contributions from n- 
phonon processes (n = 0 , 1 ,...): 


S(u, q) = S ( 0 ) (w, q) + 5 (1) (w,q) H-. (23) 


The elastic or Bragg scattering is the 0-phonon contribu¬ 
tion and is given by 

S (0) (w,q) =e~ 2W ^5{ko)NY J K^ ( 24 ) 

K 

where K is a reciprocal lattice vector, and 

e ~ 2 W ( g ) _ exp (—([ q • u(0)] 2 )t) (25) 


/cd = ( 67 r 2 ?!/) 1 / 3 ss (0.41a) 1 (21) 

is the Debye wave number. An approximate relation for 
the longitudinal mode is given by 


u;f(k) = 


e(fc) 


1 + ( k TF /k ) 2 ' 


( 22 ) 


is the Debye-Waller factor which accounts for the sup¬ 
pression of coherent scattering by thermal and quantum 
fluctuations of the ions. As mentioned earlier, the 0- 
phonon contribution does not affect electron scattering. 
At low temperature electron-ion scattering is dominated 
by the 1 -phonon contribution 


S (1) ( W ,q) 


e~ 2W 


2 M 



[q • G s (k)] 2 
w s (k) 


<5 3 (K + k-q) 


S(u> - u s ( k)) 

e w s(k)/T _ 1 


S(lo + w s (k)) 
l - e -^(k)/T » 


(26) 


where the phonon momentum k is restricted to the first Brillouin zone 13 ng. In this case S' K (q) in Eq. ([l0]) can 
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FIG. 2. (Color online) The one-phonon approximation results of S' K (q ) (thick solid curves) and S'(q) (thick dashed curves) for 
the LD matter at the three temperatures as labelled. The results of S' K (q) obtained using the fitting formula by Potekhin et al. 
[10] . which is based on the harmonic approximation for the one-component Coulomb plasma and which includes multi-phonon 
contributions, are also shown (as dot-dashed curves) for comparison. 


be replaced by 


static structure function S'(q) in OPA: 


S° FA {q) = [ dw (S (1) (w,q)\ w K (u/T,q). (27) 
J -oo ' 7 9 

Note that large-angle scattering involves a finite |K| 

|k| where the crystal absorbs a large component of the 
momentum. This is well-known as the Umklapp process 
in solid state physics [16] . Flowers and Itoh [7] realized 
that these processes dominate over normal process (with 
K = 0) in the neutron star context for typical tempera¬ 
tures of interest because pp kp>. However, transitions 
with small k and finite K are suppressed due to electronic 
band structure effects which we shall now briefly discuss. 

Although it is generally a good approximation to as¬ 
sume that electrons are free, on patches of the electron 
Fermi surface which intersect with the Brillouin zone 
boundaries, the effect of the periodic background ion po¬ 
tential is large. It distorts the Fermi surface and creates 
a band gap in the electron spectrum at the Fermi surface 
which is given by 


Ae(pp) 


4e^ b ~ w (pf) f( P f ) 

37t e(p F ) 


(28) 


S OPA (q)= r dw (5W(w,q)) • (30) 

J -oo x 7 q 


In Fig. 1 we compare S° FA (q) and S OFA (q) of the LD 
matter at three different temperatures, T/Dp = 0.03, 
0.1 and 0.3, respectively. The right panel of this figure 
shows that, even at T = 0.3 Dp, S' K (q) = S'(q) is al¬ 
ready a good approximation. The left panel of Fig. [2] 
shows that, at low temperature T = 0.03 Dp, the expo¬ 
nential factor l/(e“/ T — 1) in Eq. (151 dominates, and 
S' K (q) < S'(q) in most of the range of q. The middle 
panel of Fig. [2] with T = 0.1 Dp illustrates the com¬ 
petition between two factors in the expression of w K at 
moderate temperatures, which are discussed earlier fol¬ 
lowing Eq. |I6| ) in Section [II] For large-angle scattering 
with qa > 5 the exponential factor still dominates, and 
S' K (q) < S'(q). For small-angle scattering with qa < 5, 
however, the factor /q 2 dominates, and S' K (q ) > S'(q). 
We note that at low temperature S' K (q) has a peak at 
q —> 0 because 


where F(pp) is the charge form factor of the nucleus [T8] . 
This gap can suppress the Umklapp processes when 


T <T V ~c t Ae(p F ), 


(29) 


[1M2T] . FromEq. (29) we can deduce that Tp < 10 _2 Dp. 
In what follows we restrict our analysis to the regime 
where T is in the range 10 -2 — 1 Dp where the effects 
due to the band gap in the electron spectrum can be 
safely neglected m- 

To determine the temperature regimes where the static 
approximation S' K (q) « S'(q) is valid, we compute the 


w K (w/T < l,qa < 1) ss 1 + ^ (y) i 

and the second term in this expression always dominates 
under the typical conditions in the neutron star crust. 

For comparison we also show in Fig. [2] the results of 
S' K (q) obtained using the fitting formula by Potekhin 
et al. [ fTO] which is based on the harmonic approxima¬ 
tion for the one-component Coulomb plasma and which 
includes multi-phonon contributions. 
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IV. STATIC STRUCTURE FUNCTION AND 
MONTE CARLO SIMULATIONS 


The neutron star crust spans the regimes where purely 
classical simulations are sufficient and where quantum 
effects start to play a significant role. We can use clas¬ 
sical and Quantum Monte Carlo simulations (CMC and 
QMC) to address these conditions. The QMC simula¬ 
tions have the classical simulations as a specific limit. 
Both the CMC and QMC calculations can easily address 
the static structure function S(q). They can also be used 
to compute further information about the energy depen¬ 
dence of the response as well as other observables. 

In CMC the kinetic and potential energies are indepen¬ 
dent variables. Hence the positions of the particles can 
be sampled independently of their momentum. We use 
a simple version of the Metropolis Monte Carlo method 
to sample the positions of the nuclei in periodic bound¬ 
ary conditions at fixed density and temperature. The 
simulations use N > 1000 particles, initially at predeter¬ 
mined lattice sites in a periodic cubic box with length 
L = (N/m) 1 / 3 . Proposed particle moves {x,} —► {x'| 
have equal transition probabilities as their reverses: 


T({xJ -A {x'J) = T({x'} -A {xj), 
and they are accepted with probabilities 

e -AE/T a j£ AE > o 


P({xJ -> {x'J) = 


1 


otherwise. 


(31) 


(32) 


For CMC T a = T, and the energy change is the same as 
the change in total potential energy: 

AP = P pot ({x'})-E P ot({x i }), (33) 


fermions. Because the characteristic distance Xj of ion 
motion in lattice is much shorter than the inter-ion 
distance a, quantum statistics (the boson or fermion 
nature of nuclei) is not important, and we can consider 
the particles as distinguishable. 

In (path integral) QMC simulations the imaginary time 
or inverse temperature /3 = 1/T is split into Np slices. 
Each slice is a classical V-particle system described above 
but with effective temperature 

= (At) "‘ = (i) 1 ■ (35) 

For QMC each imaginary-time slice involves a high- 
temperature expansion of the propagator exp(— HAt). 
For a large enough number of slices Np the results are 
independent of the number of slices. Typically of order 
10 slices are required in the present calculations. 

As in CMC, the Markov chain is again constructed by 
moving the particles according to the acceptance proba¬ 
bility defined in Eq. ( |32| . For QMC the energy change 
includes the changes in both kinetic and potential ener¬ 
gies: 


A E — [E l pot({x,; (T }) + -Ekin({Xj )Cr })] 

[Ppot ({ x i,cr } ) 4“ 7^kin({ x i,a'})] i 


where 


N Np 


Pkin({^i,cr }) — EE 


( X 'i,CT + l X i,o-)° 

2M (At) 2 

Np 

ppot({ x ij <7}) = y (y ( v(\Xi t(J — x y, cr 1) 

a=l i<j 


(36) 


(37) 

(38) 


where 


£pot({ x J) = E y (l Xi_x jD- ( 34 ) 

i<j 

In Eq. ( |34[ ) the sums over i and j run over the parti¬ 
cles in the simulation volume plus their periodic images. 
Typically plus or minus one image in each direction (i.e. 
27 periodic boxes in total) is sufficient in these simu¬ 
lations because of the screening of the ion-ion potential. 
Standard Ewald summation is also possible but would be 
slower. Detailed balance ensures that the Markov chain 
constructed with the method described above will con¬ 
verge eventually to sample particle positions proportional 
to the partition function. 

Quantum fluctuations become important when 
T/flp < 1. For such scenarios we used path integral 
QMC simulations (see, e.g., [22]). The single position 
for each particle in the classical simulation becomes a 
path in path integral simulations with periodic bound¬ 
ary conditions in imaginary time. Boson or fermion 
path integrals would require us to include exchanges 
in the imaginary time boundary conditions with the 
appropriate statistics, e.g. —1 for odd permutations of 


with ^i-i.Np+i — x i.i- Clearly, a CMC simulation can 
be considered as a special case of QMC simulation with 
Np = 1. 

The static structure function is then obtained from the 
points sampled after convergence. In Monte Carlo simu¬ 
lations 


S(q) = 


NN R 


I Np N 

e iq 

\<7—1 ij = 1 


(39) 


q ,t 


which includes both the one-phonon and multi-phonon 
contributions. Because of the periodic condition, 


27T „ „ 

q = -y-(n x *L + n y y + n z z) 


(40) 


take discrete values, where n x ^ ViZ ') are integers. To ob¬ 
tain the inelastic part of the static structure function 
S'(q) we simply remove the points that correspond to 
the Bragg peaks in the BCC lattice. Other detailed struc¬ 
tures predicted by QMC and CMC, i.e. the smaller peaks 
and troughs away from the Bragg peaks (see Fig. [ 3 ]), are 
finite-size artifacts whose amplitude decreases with in¬ 
creasing particle number in the simulation. However, the 
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FIG. 3. (Color online) The inelastic part of static structure function S'(q) of the LD matter at the two temperatures as labelled. 
The results are obtained using the one-phonon approximation (thick dashed curves), and quantum and classical Monte Carlo 
simulations (solid and dotted curves), respectively. 


integrated strength over any reasonable interval in q is 
physically relevant and is insensitive to finite-size effects 
after the numerical convergence has been achieved. 

At high temperature T > Dp all phonon modes are ex¬ 
cited, S CMC (q) and S^ MG (q), which are S'(q) obtained 
using CMC and QMC simulations, respectively, should 
agree. But at low temperature T -C Dp quantum fluc¬ 
tuations become prominent, and S^ MG (q) > S GMG (q). 
In Fig. [ 3 ] we compare S GMG (q) and S ,< ^ MC (g) for the LD 
matter at two different temperatures, T/Dp = 0.1 and 
0.56, respectively. Indeed, S GiMG (q) is clearly larger than 
S GMG (q) at T = 0.1 Dp, but it is somewhat surprising 
that S GMG (q) and S ,QMC (( 7 ) agree very well even at tem¬ 
perature as low as T = 0.56 Dp. 

For comparison we also show S OPA (q) in Fig. [3 This 
figure shows that S GiMC {q) and S OPA (q) agree with each 
other at low T and/or small q, although there exist rapid 
oscillations in Monte Carlo results because of the finite 
size of the system. At high T and/or large q multi¬ 
phonon contributions are signfiant, and the one-phonon 
approximation breaks down. 

Note that in Eq. (39) only the equal-time correlator has 
been evaluated. By including an offset in the imaginary 
times between the evaluations of positions of particles i 
and j, 


S(sAr,q) 


1 

NNp 



e iq(xi, CT -Xj,„ +s ) 



, ( 41 ) 


one can obtain information about the energy dependence 
of the response [22j. It is also possible to calculate the 
properties of MCP in both the classical and quantum 
regimes. This would require simulations significantly 
larger than the OCP studied here, to ensure that the 
periodic boundary conditions do not impact the results. 
Simulations of this magnitude should be readily achiev¬ 
able on modern parallel computers. 


V. RESULTS AND DISCUSSION 


We have calculated the thermal conductivity of OCP 
for the LD and HD ambient conditions outlined in Ta¬ 
ble |T] for the temperatures of interest to neutron star phe¬ 
nomenology. Our primary motivation to study the simple 
one component system was to obtain a quantitative un¬ 
derstanding of the effects of dynamical information and 
quantum fluctuations at intermediate temperatures in 
the range 0.1 Dp < T < Dp. To this end we use the 
various approximate methods outlined in previous sec¬ 
tions to calculate S' K (q), which is the kernel function for 
computing effective electron collision rate v K in Eq. (10). 
We then calculated n for the catalyzed neutron star mat¬ 
ter with densities 10 1Q g cm -3 (LD) and 10 11 g cm -3 
(HD), respectively. The results are shown in Fig.[4]where 
the thermal conductivity is obtained by replacing S' K (q) 
with S GPA (q) (thick solid curves), S OPA (q) (thick dashed 
curves), S GPA (q) with simple phonon dispersion relations 
[see Eqs. ( |20| ) and ( [22] )] (thin solid curves), the fitting 
formula of S K (q) based on the harmonic approximation 
[TO] (thin dot-dashed curves), S^ MG (q) (filled circles) and 
S GMG (q) (filled squares), respectively. 

A careful comparison of the results obtained using dif¬ 
ferent approximations for the function S' K (q) provides the 
following useful insights: 


1. We find that it is adequate to set S' K (q) = S'(q) in 
calculating n at temperature as low as T ss 0.1 Dp. 
A comparison between the thick dashed curves ob¬ 
tained using S OPA (q) and the thick solid curves ob¬ 
tain using S GFA (q) supports this conclusion. The 
validity of this approximation is expected and well 
known for T > Dp. Our results show that S' K (q) ~ 
S'(q) even at T = 0.3 Dp (see Fig. [2j). Our results 
also show that, for 0.1 < T /Dp < 0.3, this approx¬ 
imate method can still be used to compute k even 
though S' K (q) and S'(q) differ. This is because there 
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FIG. 4. (Color online) Thermal conductivity of the LD (left panel) and HD matter (right panel) in units of 10 18 erg cm _1 s~ 1 K _1 
and as a function of temperature. The results are obtained by replacing S' K (q) in Eq. ( |10[ ) with S° PA (q) (thick solid curves), 
S OPA (q) (thick dashed curves), S® PA {q) with approximate phonon dispersion relations {see Eqs. ( |20[ ) and ( |22[)] (thin solid 
curves), fitting formula of 5")(q) for one-component Coulomb plasma based on the harmonic approximation [10] (thin dot- 
dashed curves), 5 ,CJMC (q) (filled circles) and , 5 CMC (q,) (filled squares), respectively. 


are two competing factors in w k (lo/T , q) which are 
discussed earlier following Eq. (30) in section |TIJ 


2. Multi-phonon effects become relevant for T /T m < 4 
in our simulations. At higher T or lower tempera¬ 
ture the one-phonon approximation is adequate for 
OCP but is sensitive to the phonon dispersion re¬ 
lation. This is evident when we compare the thick 
solid curves, the thin solid curves and thin dot- 
dashed curves, which are obtained using the ex¬ 
act phonon dispersion relations from the dynami¬ 
cal matrix, the approximate phonon dispersion re¬ 
lations, and the harmonic approximation method 
with multi-phonon contributions m, respectively. 


3. The comparison between the results obtained using 
CMC and QMC simulations, shown by the filled 
squares and circles, respectively, indicates that 
quantum effects in thermal conductivity are signif¬ 
icant when T < 0.3 Up where classical calculations 
systematically underestimate S'(q). At T « 0.1 flp 
CMC results overestimate n by 40 — 50%. 


4. The fitting formula for S’ K {q) based on harmonic 
approximation [TP] (dot -dashed curves) works quite 
well for the one-component Coulomb lattice in the 
whole temperature range which we have studied. 
This can be seen when we compare them with the 
thick solid curves at T < 0.3 flp and filled cir¬ 
cles/squares at T > 0.3 flp which are obtained 
using the one-phonon approximation and Monte 
Carlo simulations, respectively. At the highest tem¬ 
perature where T « T m the results obtained using 
the harmonic approximation include multi-phonon 
excitations and agree well with the QMC results. 
This indicates that anharmonic effects are small 


even in this regime. At the lowest temperatures, 
although the S' K (q ) obtained from the fitting for¬ 
mula based on the harmonic approximation differs 
from that obtained in the OPA (see Fig. [2|, the 
predictions for the thermal conductivity agree well 
as already discussed in nnj. 

Some of the trends emerging from these comparisons 
could have been be expected qualitatively. As previously 
alluded to, this systematic quantitative comparisons be¬ 
tween QMC results and those obtained using the stan¬ 
dard electron-phonon treatment provide a basis to asses 
the viability of using QMC calculations of S(q) at low 
temperature for complex multi-component systems. In 
the standard treatment, multi-component systems are 
modeled as a regular lattice plus uncorrelated impurities, 
and electron scattering is assumed to arise due to inco¬ 
herent contributions from electron-phonon and electron- 
impurity scattering. This treatment fails when the spa¬ 
tial distribution of the minority species is correlated, and 
QMC is a viable technique to calculate the role of these 
correlations in strongly coupled plasmas with T 1 in 
the regime when T < flp. 

Our results show that S(q) obtained from QMC is ad¬ 
equate to calculate k for T > 0.1 flp, and CMC may 
be adequate to compute thermal conductivity of OCP at 
T > 0.3 flp. For lower temperatures, more detailed in¬ 
formation about the energy dependence of the response 
is needed and we have briefly commented on how we can 
accesses this in the discussion following Eq. f4l] ) . With a 
modest increase in computing resources, we have found 
that QMC can be used to study large multi-component 
plasmas, and we are currently pursuing calculations of k 
for ambient conditions encountered in accreting neutron 
stars. These results will be reported elsewhere. 
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